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INTRODUCTION 


During  the  past  few  years  many  research  groups  have  engaged  in 
study  of  the  performance  characteristics  of  a  plug  nozzle.  As  yet  a  com¬ 
puter  program  to  study  the  flow  pattern  and  performance  has  not  been 
reported.  This  report  summarizes  a  basic  analytical  method  and  describes 
a  computer  program  based  on  this  method. 

The  basic  characteristic  equations  were  derived  by  assuming 
rotational  flow,  so  that,  in  the  future,  shock  equations  could  be  added  to 
the  present  calculations  without  difficulty.  The  gas  is  assumed  to  be 
perfect  and  inviscid.  Friction  loss  on  the  nozzle  wall  is  ignored,  and  the 
base  pressure  of  the  plug  is  computed  by  using  Korst's  theory. 

The  present  numerical  method  has  been  programmed  in  IBM  7040 
FORTRAN  IV,  and  two  sample  calculations  are  presented  in  this  report. 
This  program  can  be  used  to  examine  the  performance  of  various  plug 
nozzle  design  concepts. 
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ANALYSIS 


The  flow  field  of  a  plug  nozzle  is  formed  by  an  axisymmetric 
internal  plug  with  an  external  solid  boundary  at  the  upstream  and  free 
expansion  at  the  downstream.  It  consists  of  a  base  pressure  region  at 
the  end  of  the  plug  if  the  plug  is  truncated. 

The  method  of  characteristics  is  used  to  calculate  the  supersonic 
flow  fields  and  the  Prandtl-Meyer  relations  are  used  to  calculate  the  flow 
properties  of  the  lip  of  shroud.  The  base  pressure  problem  is  soived  b y 
using  Korst's  theory.  The  gas  is  assumed  to  be  perfect,  inviseid,  and  the 
flow  field  is  assumed  to  be  steady,  rotational  and  axisymmetric. 

Basic  Equations  of  the  Method  of  Characteristics 

The  characteristic  equations  for  axisymmetric,  steady  and 
rotational  flow  vised  in  this  analysis  were  presented  by  A.  II.  Shapiro  in 
Reference  1.  The  characteristic  equations  were  derived  from  continuity, 
energy  and  Euler's  equations.  The  detailed  derivations  were  also  shown 
in  Reference  2.  There  are  two  families  of  characteristics: 


The  geometric  properties  of  the  characteristics  provide  other  relations: 


Left  Running  Characteristic 


dy 

dx 


=  tan  (0  +  (3) 


(3) 


Right  Running  Characteristic 


dy 

dx 


=  tan  ( G  -  (3) 


(4) 


Writing  Equations  3  and  4  in  finite  difference  form  and  solving  for  x,  y, 
one  obtains: 


The  last  terms  in  Equations  1  and  Z  are  to  take  into  account  the 
entropy  change  in  the  flow  field.  In  order  to  compute  the  entropy  change 
along  a  characteristic,  the  entropy  is  assumed  to  be  constant  along  a 
streamline  and  varied  across  a  streamline.  Since  the  entropy  gradient 
is  not  large,  it  is  also  assumed  to  be  constant  in  each  small  region. 

The  derivations  were  presented  in  Reference  1  and  the  expression 
can  be  written  as  follows: 
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{ S  2  -  Si)  (x3  -  Xj) 


sin  p 

cos  (0+| 3) 


s3  =  s,  + 


(x3  -  Xj) 


sin  P 


+  (x3  -  x2) 


cos(0  +  /3)J  3"  Icos(0-j3) 


sin  P 


The  velocity  and  the  flow  angle  at  point  3  can  be  solved  by 
combining  Equations  1  and  2: 


V,  - 


cot  P  cot  P 

~~V~  V~ 

-J13  C  _l23 


e2.o1+(^i£)  v1+  (HL£) 


sin  p  sin  0 
y  cos  (0  +  p) 


,  .  sin  p  sin  0 

(x3-x,)+  - - —  X. 

L  y cos  (°  -  V)  J,3 


-  X3) 


T  T 

—  sin  p  cos  p  (s3  -  s  i)  -  —  sin  p  cos  p  (si  -  s2) 

La2  J 13  La2  J-3 


03  =  0i  +  (V3  -  V,) 


cot  p  sin  p  sin  0 

V  _  y  cos  ( 0  +  p) 


(X3  -  X|) 


+  —  sin  p  cos  (s3  -  s 

_  a2  J  i3 


When  a  right  characteristic  intersects  the  boundary,  as  shown  in 
Figure  2,  the  intersection  can  be  solved  by  the  following  equations: 


yB  =  yBl  +  (xB  "  XB1)  tan  G1 


VB  =  Yi  +  (XB  -  xi)  ttan  (0  -  /^)J  j 
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Figure  1.  Nomenclature  for  Method  of  Characteristics  in 
Rotational  Flow  Field  Calculations 


Figure  Z.  Nomenclature  for  Method  ol  Characteristics  in 
Boundary  Point  Calculation 
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where 


In  order  to  compute  the  flow  properties  at  the  end  point  of  the 
boundary,  it  is  necessary  to  insert  a  characteristic  at  that  point  as 
shown  in  Figure  3.  When  a  left  running  characteristic  is  inserted,  the 
intersection  may  be  solved  by  using  the  following  equations: 


,  y<i  -  y3 

[tan  (0  +  /3 ) ]  = -  (17) 

34  x4  -  x3 


y 4  -  y i  _  y2  -  yi 

X  4  “  X  |  XC  ^  Xk.  | 


(18) 


Solving  x4  from  Equations  17  and  18  yields 


x4  = 


yi  -  x3  +  x 3  [tan  (0  +  |j)]  -  x,  (— - —  ) 

44  \  x  i  -  x  J  / 


[tan  (0  +  /])] 


34 


Vi  -  yi 

X2  -  X, 


(19) 


The  veloc  ity  at  point  3  can  lie  computed  by  using  the  following  equation: 

sin  (1  sin  0 


V3  =  V4  +  [V  tan  jj] 


34 


(03  -  04)  + 


L  y  cos  +  ]S) 


(x3  -  x4) 


34 


T  .  „ 

(  S  3  -  S  4 ) 

—  sm  p  cos  p 

,  a2 

34 

(20) 


Similarly,  the  relation  for  a  right  running  inserted  characteristic  can  lie- 
written  as  follows: 


yi  -  y-s 


+  x 3  [tan  (0  -  /i)l  -  X]  - — ) 

44  \  X  i  -  X  ]  ' 


X  4  = 


[tan  (0  -  /$)] 


34 


Yz  -  yi 

X2  -  X] 


(21  ) 
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(x3  -  x4) 


V3  =  V4  +  [V  tan  0] 


34 


r 

(e4  -  e3>  + 

sin  0  sin  0 

_  y  cos  (0  -  0)_ 

34 


T 

—  sin  0  cos  0 

(s  3  -  s4) 

_  a2 

34 

Transonic  Region 

The  transonic  flow  near  the  throat  of  a  nozzle  requires  special 
treatment  because  the  method  of  characteristics  is  not  valid  in  this 
region.  The  Sauer  analysis  in  Reference  3  offers  a  solution  to  this 
problem.  The  solution  was  presented  as  a  power  series  and  the  derivation 
was  based  on  two  dimensional,  small  perturbations  theory. 


u  =  ox  -t 


Y+l 


z  —  i  , 
o  y  + 


(Z2) 


(Y+l 


o  x  y 


( v  n  V 


y  ‘ 


(2  i) 


where 


O' 


:  \  (y+D 


Ps  Vs 


U4) 


The  values  ps  and  ys  can  be  obtained  from  the  geometry  of  a  nozzle  as 
shown  in  Figure  4. 


Base  Pressure  Region 

When  a  plug  nozzle  is  truncated,  the  base  pressure  becomes  an 
important  parameter  affecting  the  nozzle  performance.  Korst's  analysis 
in  Reference  4  provides  an  approach  to  this  problem.  The  derivations  are 
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based  on  two  dimensional  turbulent  flow  with  constant  pressure  mixing. 
The  essential  feature  of  the  flow  model  is  shown  in  Figure  5.  The 
boundary  layer  at  separation  is  assumed  to  be  thin  compared  to  the  length 
of  the  jet  mixing  region  and  no  n.ass  is  assumed  to  bleed  into  the  wake. 
The  Crocco  number  is  defined  as  follows: 


JL  + 

y!  1 


For  isoenergetic ,  fully-developed,  turbulent,  constant  pressure, 
jet  mixing  profiles,  the  velocity  ratio  is 


5  =  <* +  urf  ,l) 


where 


erf  t|  = 


j  c-H'd/* 


T]  -  O’  •- 
X 


<r  =  1  2  +  2.  758  M2 


In  the  case  of  no-bleed  into  the  wake,  the  Crocco  number  at  j  streamline 


cd'  =  ^  c/ 
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Assuming  M2a  -  the  isentropic  relation  gives 


P 4 

P3 


1 


(1  - 


Cd2^-1 


(31) 


The  flow  turning  angle  ]02  can  be  obtained  from  the  Prandtl  -Meyer 
relation  by  assuming  jG2  ~  3O4  ,  and  the  Mac!)  number  at  region  1  can 
also  be  obtained  from  the  Prandtl  Meyer  relation: 


M  ,  -  f  (v  ,)  (33) 

where 

v  1  v  i  -  |02  (33) 


Using  isentropic  relations,  the  base  pressure  for  back  step  can  be 
computed  as  follows: 


Pi 

Pb 


M  , 


(  M) 


By  assuming  a  series  of  values  for  M2  in  Equation  25,  and  carrying 

pb 

through  the  whole  procedure,  a  curve,  vs  Mj  can  be  obtained. 

P 1 

In  order  to  take  into  account  the  effect  of  boattailing,  Korst's 
Reduced  Mach  Number  Concept  has  to  be  used  to  extend  the  previous 
technique.  0a''  was  defined  as  a  streamline  angle  at  which  M  =  1  pro¬ 
duced  by  the  Prandtl  -  Meye  r  relation  from  M  ia  and  0  ja  . 


M  ia  -  M|a  (-0ta  -  0a  ) 


1  1 


(35) 


Figure  5.  Korst's  Flow  Model 


Figure  6.  Geometrical  Configuration  of  Base  Region 
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A  reduced  Mach  number  can  be  determined  as  follows: 


^■red  ~  ^red  ®a  ) 


(36) 


Pb  p2 

— —  vs  Mj  curve  obtained  from  the  previous  technique  is  taken  as  — - 

pi  pred 


vs  Mrccj.  Then,  the  base  pressure  can  be  computed  by  using  the 
following  relations: 


pb  _  P2  poa  Pred 

Pi  Pred  Pi  Poa 


(37) 


where 


P) 

P  oa 


Pi 

Poa 


(M,) 


(38) 


Pred 

P  oa 


Pred 


oa 


(Mred) 


(39) 


Numerical  Procedure 

The  computations  consist  of  several  distinct  parts:  the  calcu¬ 
lations  of  a  starting  line,  field  points  and  boundary  points,  Prandtl- 
Meyer  expansion,  and  the  base  pressure.  The  starting  line  is  determined 
by  using  Equations  22,  23  and  24.  The  computation  of  field  points  and 
boundary  points  is  performed  by  a  regular  iteration  scheme.  The 
coefficients  of  mean  values  are  employed  in  the  process  as  suggested  by 
Darwell  in  Reference  5.  The  calculation  of  Prandtl  -  Meyer  expansion 
takes  part  in  the  process  when  the  last  upper  boundary  point  is  obtained. 
When  the  last  point  of  the  lower  boundary  is  reached,  the  bust:  pressure 
Computation  is  employed. 
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The  cumulative  vacuum  thrust  is  made  up  of  the  momentum  flux 
and  the  pressure  thrust  at  the  starting  line  plus  the  pressure  integral  on 
the  boundaries.  The  mass  flow  rate  across  the  segment  12  as  shown  in 
Figure  7  is 


m  =  p  12  V,,  A  12  cos  (4>  -  012) 


where 


4>  =  tan" 


X]  -  x2 

ya  -  yi 


A  ]2  =  TT  (y  1  +  y2)  ^{x2  -  X))2  +  (yi  -  yz)‘ 


The  momentum  flux  and  pressure  thrust  at  the  segment  12  at  the  starting 
line  are 


M0  =  —  V12cos  0|2i  F|>  A^  cos  4 
g 


The  pressure  integral  at  segment  12  on  the  plug  is 


Mp  =  P  12  A  12  cos  4> 


The  cumulative  vacuum  thrust  can  be  computed  as  follows: 


M0  + 


The  vacuum  thrust  coefficient  is  defined  as  follows: 


(Cp) 


T  +  ”  r] 
PQ  Arp 


The  numerical  procedure  described  in  this  report  has  been 


programmed  in  IBM  7040  computer  FORTRAN  IV  language. 
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REMARKS  ON  CALCULATIONS 


The  accuracy  of  the  present  method  depends  on  the  net  size 
chosen  for  the  calculations.  In  other  words,  the  smaller  the  net  size 
one  chooses,  the  more  accurate  the  results  one  can  obtain.  When  the 
small  net  size  is  used,  of  course,  the  points  used  to  describe  the  contour 
should  be  more  accurate.  There  are  two  ways  to  control  the  net  size. 

One  is  to  control  the  number  of  points  at  the  starting  line,  and  the  other 
is  to  control  the  number  of  rays  at  the  lip  of  a  nozzle. 

When  the  inclined  angle  of  the  lower  wall  becomes  large,  Un¬ 
reduced  Mach  number  computed  froin  Equation  36  differs  from  the  Mach 
number  at  the  edge  in  a  great  amount.  Tins  difference  may  cause  the 
base  pressure  to  be  greater  than  the  pressure  on  the  boattailed  portion  as 
shown  in  Equation  37.  In  this  case,  it  may  indicate  separation  and  the 
theory  becomes  invalid. 
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SAMPLE  RESULTS  AND  DISCUSSION 


The  program  has  been  used  to  compute  several  test  cases.  Two 
typical  cases  are  selected  for  presentation  in  this  report.  An  external 
expansion  plug  nozzle  was  designed  by  using  the  program  in  Reference  6. 
In  order  to  compute  a  starting  line  for  the  analysis,  the  simple  wave 
relation  was  employed.  The  computer  results  are  shown  in  Figure  8. 

The  vacuum  thrust  coefficient  is  about  one  percent  higher  than  the  design 
value,  but  the  design  method  was  assumed  as  a  simple  wave  throughout 
the  whole  flow  field.  An  inter  na  1  -  externa  1  expansion  plug  nozzle  was  also 
computed.  The  result  and  the  flow  pattern  are  shown  in  Figure  9. 

When  a  nozzle  contour  is  not  well  described  or  a  compression 
region  occurs  in  the  flow  field,  the  characteristics  overlap  indicating 
that  a  shock  is  being  formed  in  that  region.  Il  the  shock  is  weak,  the 
present  program  carries  on  the  calculations  by  assuming  an  isentropic 
process.  A  shock  routine  must  be  developed  to  analyze  a  nozzle  with  a 
strong  shock.  The  Rankine - Hugoniot  equations  are  normally  used  for 
this  purpose. 

In  the  derivation  of  the  transonic  solution,  the  second  degree  of 
the  velocity  components  was  ignored.  Therefore  a  significant  amount  of 
error  would  be  introduced  to  the  result  if  the  Mach  number  of  the  starting 
line  were  high.  In  the  case  shown  in  Figure  9,  two  percent  of  error  in 
vacuum  thrust  was  found  when  the  initial  Mach  number  changed  from 
1  .  05  to  1.15. 
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Figure  8.  Flow  Field  of  an  External  Expansion  Plug  Nozzle 
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Figure  9.  Flow  Field  of  an  Internal-External  Expansion  Plug  Nozzle 


This  program  is  suitable  for  a  basic  study  of  plug  nozzle 


performance.  In  order  to  improve  the  quality  of  the  result,  the 
following  items  are  recommended  for  future  work. 

1.  To  develop  a  shock  routine  there  will  be  no  difficulty  because 
rotational  flow  was  assumed  in  the  present  program. 

2.  To  include  real  gas  equations  in  the  computation. 

3.  To  take  into  account  the  friction  loss  on  the  nozzle  walls. 
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DESCRIPTION  OF  DATA  INPUT 
AND  OUTPUT 


Input 

This  program  requires  the  following  input  data: 


(1)  Nozzle  components 


FE  (4)*) 


ROS  (ps) 


YS  (ys) 


GAM  (Y) 
XM  (Mest) 
P  (P0) 

T  (T0) 

RT 

R 

N 

N1 

N2 


throat  plane  inclined  angle;  degrees  for  internal- 

external  expansion 

radians  for  external  expansion 

used  only  for  internal-external  expansion 
equals  0.  for  external  expansion 

radius  of  nozzle  throat;  used  only  for  internal- 

external  expansion 

equals  0.  for  external  expansion 

ratio  of  specific  heats 

initial  Mach  number 

total  pressure 

total  temperature 

radius  from  the  throat  to  origin,  used  only  for 
interna  1- external  expansion 
equals  0.  for  external  expansion 

gas  constant 

number  of  points  on  starting  line,  must  be  <  100. 
number  of  lower  wall  contour  points,  must  be  <100. 
number  of  upper  wall  contour  points,  must  be  <100. 


(2)  A  title  or  job-description  card 

(3)  N1  lower  wall  contour  points  given  as  Cartesian  coordinates 


(4)  N2  upper  wall  contour  points  given  as  Cartesian  coordinates 


(5)  KK  --  If  input  is  in  feet  kk  =  0 

--  If  input  is  in  inches  kk  =  1 

(6)  KODE 

(a)  If  KODE  =  1,  read  starting  line  for  an  internal-external 
plug  nozzle  expansion 

(b )  If  KODE  =  2,  compute  starting  line  for  an  internal- external 
plug  nozzle  expansion 

(c)  If  KODE  =  3,  compute  starting  line  for  an  external  plug 
nozzle  expansion 

(7)  KODE  used  for  external  expansion  only 

(a)  If  KODE  =  1  use  standard  starting  line  calculations 

(b)  If  KODE  =  2  use  a  special  option  in  calculating  the  starting 
line:  Mest  =  Mp  and  -  Op. 

(8)  NU  (r|)  the  number  of  corner  rays  to  be  computed,  must  be 

<100. 

PA  (Pa)  --  Ambient  Pressure 
KKD  --  If  KKD  =  0,  PA  is  in  lbs/sq  ft 
--  If  KKD  =  1,  PA  is  in  lbs/sq  in 

Input  cases  can  be  stacked  and  processed  several  at  a  time.  If  a 

bad  data  case  is  found,  the  remaining  data  cases  will  not  be  processed. 

This  is  due  to  the  computer  system,  not  the  program. 


INPUT  DATA  AND  FORMAT 


1 

CO 

:  <b 

•M 

u 

03 

+-» 

4-> 

u 

P 

p 

03 

a 

cx 

c 

c 

u 

►H 

Wh 

u 

c 

O 

o 

•  H 

0 

*— * 

- — - 

(l) 

•H 

r— ( 

t—i 

c 

c 

CX 

p 

•  H 

c 

c 

£ 

•  H 

•H 

03 

u 

TD 

"0 

rC 

U) 

(U 

<D 

a 

<D 

G 

c 

rH 

Q 

•rl 

•  »H 

03 

0) 

<u 

00 

TJ 

TD 

r- 

<b 

<U 

Sh 

^  £ 

U  C 

0  _ 

03  0 

rt  0 

c  00 

•r- < 

•  »H 

fc  t*- 

«  *; 

— «  a. 

^  a. 

!  P  ' 

O  -H 

o  t1 

r*  r— < 

r 

,£> 

X)  * 

•  fH 

("  o 

r  o 

X  r— < 

C  in 

C  in 

03  O 

>,  ai 

>■  <u 

to  T3 

C/3  ■v 

ZS 


PA,  KKD  I  Symbols  are  defined  in  (8)  of  Input 


OOOE-f-Ol  +  O.  OOOOOOOOE+OO-  0.  00000000E+00+0.  1  4000000E+0I +0.  12000000E+01 


cc 


Output 


(1)  Units  of  variables 

(2)  Job  title 

(3)  Input  conditions 

(4)  Upper  wall  contour 

(5)  Lower  wall  contour 

(6)  Starting  line  points 

X  Y  M  THETA  T  P 

*  1  1  1  11 

where  X,  Y  are  Cartesian  coordinates;  M  is  Mach  number, 
THETA  is  flow  angle,  T  is  temperature  (°R),  and  P  is 
pressure 

(7)  Internal  expansion 

(a)  Field  routine  points 

X  Y  M  THETA  T  P  ITR 

111  1  111 

where  ITR  is  the  number  of  iterations  before  convergence  in 
calculations 

(b)  Body  point  routine  point 

X  Y  M  THETA  T  P  ITR 

Field  and  body  points  alternate  until  the  last  point  on  the 
upper  wall  contour  is  reached 

(8)  External  expansion 

(a)  Insert  point 

X  Y  M  THETA  T  P 

(b)  Corner  point 

X  Y  M  THETA  T  P 
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(c)  Right  running  characteristics 

X  Y  M  THETA  T  P 

ill  1  1  1 

(d)  Field  routine  points 

X  Y  M  THETA  T  P  ITR 

111  1  ill 

(e)  Body  point  routine  point 

X  Y  M  THETA  T  P  ITR 

Field  and  body  points  alternate  until  the  last  point  on  the 
lower  wall  contour  is  reached  or  until  the  network  is 
completed. 

(f)  Insert  point 

X  Y  M  THETA  T  P 

(g)  Corner  point 

X  Y  M  THETA  T  P 

(9)  Thrust  distribution  along  the  plug 

(a)  SUMM 

(b)  CFI 

(c)  Mass  flow  rate 

(d)  X  Y  T  CF 

i  i  i  ; 

for  each  point  on  the  lower  wall 

(e)  AT  --  throat  area 

(f)  PB  --  base  pressure 

(g)  TVAC  --  vacuum  thrust 
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(h)  CF VAC  --  Vacuum  thrust  coefficient 

(i)  THRUST  --  real  thrust 

(j)  CF  REAL  --  real  thrust  coefficient 

(k)  End  of  job 
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LIST  OF  FORTRAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  BY  USING  THE  METHOD  OF  CHARACTERISTICS 


MAIN  PROGRAM 

DIMENSION  YPI400) ,  XPI400) , THI400) ,XMPt400) , TPI400) , PPI400) , 

IRXMI 200) ,RTH( 200) ,RTP(200) ,RPP(200) ,VLP(400) 

DIMENSION  XB1 (  100)  ,  X82 (100),YB1(100),YB2{10U) 

DIMENSION  FRX ( 50) , FRY  I  50) , FRV { 50) ,FRT l 50) , FRP ( 50) , FRTH{ 50) , 

I FUX ( 100),FUY(100),FUP(100),FLX{200),FLYt200) , FLP( 200) 

DIMENSION  ZC ( 100 ) , ZJ{ 100) ,XM1 ( 100) ,PBP1  ( 100 ) 

COMMON  YP,XP, TH,XMP,TP,PP,RXM,RTH,RTP,RPP,VLP,RGS, YS,GAM,GM1,G, 
1XM,N,P,T,K,L,M, J,N2,XXB2,YYB2,NU,KNT,GPl,FL,RT 
1  ,PA 

COMMON ZC,ZJ,XMl,P8 PI, NO 
RE  AD ( 5 i 52 ) NG 

READ! 5 ,1003) ( Z J  II  )  , ZC ( I ) , I  =  l , NU ) 

78  READ! 5, 100 1 ) FE , ROS , YS , GAM , XM , P , T , RT , R , N , N 1 , N2 
300  FORMAT ( 1H1 ,54X ,22HPLUG  NOZZLE  AN AL YS I S / 1  HO , 44X , 4 3HBY  USING  THE 
l  METHOD  OF  CHARACTEKISTICS////1H0, 1  OX , 5HUN I T S/ // 1H0 , 10X, 16HC00R 
1 D I  NATE  S  X,Y» 14X ,4H IN.  / 1H0 , 10X , 25H I NCL I  NED  THROAT  ANGL E ( F E J , 5X , 7H 
10EGREES/1H0, 10 X, 8HPRES SUR E , 22 X , 9HL BE / IN* IN/ 1H0, l OX, 11HTEMPERATURE, 
1  19 

IX, 1 5HDEGREE S  R ANK I NE / 1 HO , 1 OX , 1 6HG A S  CONSTANT  I R ) , 1 4 X , 26HFT  LBF/LBM 
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LIST  OF  FORTRAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  BY  USING  THE  METHOD  OF  CHARACTERISTICS 


l  DEGREES  RANKINE/ I HO , I OX , A HARE A , 26X , 5H I N* I N/ I HO ,  1 OX , 6HTHRUS T , 2 A  X , 3 
IHLBF  J 

3CI  FORMAT { I 3A6 ) 

302  FORMAT  I  1H0, I3A6) 

3C3  FORMAT  [  lFiO,  I  OX,  I  7H INPUT  CONDI  Tl CnS/ / / 1  HO , 10X , 3HF E  =  E 1 5 . 8/ 1H0,  10X, 3 
1  HR  T  =  E l 5 . d/ 1  HO ,  l  OX,  3HYS  =  EI5.8/  I  HO  ,  1 0  X  ,  6  H  R  H  0  S  =  E  1  5.8/ l  HO,  10X,6HGAMMA= 
IE  I 5. 8/1 HO, I OX, 5HMEST=E15.8/ IHO ,  LOX , 2HR  =  C I  5 . 8/ I  HU , l 0 X , 3HP0= E 15 . B / 1 H 
10,  10X, 3HT0=L15.B) 

READ(.5,30l  )  A  I  »A2»A3»A4,A5»A6,A7  ,AK,a9,A10,A1I,A12,AL3 
WR  l  I  E ( 6 , 300 ) 

WRI  TE( 6, 3U2J A  1,A2, A3, AA,A5,A6, A/, AB,A9,A10,A1 1 ,A12,A13 
READ! 5 , 1003 ) t  X81 ( I  )  , YBI II ) , I  =  I ,N1 ) , ( XB2(I) , YB2 (I) , 1=1 ,N2) 

WRIT E(C, 303)Fl,RT,YS,ROS,GAK,XM,K,P,T 
WRI  10(6,1006) 

WR  1  1  E ( 6 , 1007 )  (  X82  (  I  ) , YB2  l  I ) , I  =  1 , N2 ) 

WRITE(6,100H) 

WR  IT t ( 6, 100  7 )  ( XB1 U )  , YBl (1) , 1  =  1 ,N1 ) 

READ ( 6 , 62 ) KK 
IFIKK.EO.OJGO  TO  401 

P=P*144. 
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LIST  OF  FORTRAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  BY  USING  THE  METHOD  OF 


YS=YS/12. 

ROS=ROS/ 1 2 . 

R1 =RT/  12 . 

DO  A 02  K  =  1  *  M 
XB1IK)=XB1 ( K ) / 12 . 

4C2  YB  1  (  K  )  =  Yb  L  IM  /  1 2  . 

DO  403  K=1,N2 
XB2 I K ) -XB2 { K  )  /  12. 

403  YB^IKJ=YB2IK)/12. 

4 C  l  XXts2  =  Xu2IN2) 

YYU2=YB2(N2) 

10C6  FORMAT { LHO, 18HUPPER  WALL  CONTOUR/ 1H0 , 7X , 1HX  ,  16X,  1HY ) 
10C7  FORMAT ( IhC ,2 ( E IS. 8 ,2X ) ) 

1 0C8  FORMAT ( 1H0 , lbHLOWER  WALL  CONTOUR/ IHO , 7 X , LHX , 1 6 X , 1 HY ) 
NF  =  2 
NU  =  0 
KN  T  =  l 

GM 1 =OAM- 1 . 

GP 1 =G AM+ 1 . 

G  =  32 . 2 


CHARACTERISTICS 
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LIST  OF  FORTRAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  BY  USING  THfc  METHOD  OF  CHARACTERISTICS 


WR 1 rt ( 6, 1002 ) 

R£AD(5,52)KODb 

C  K0UE=1 — REAu  START  LINt 

C  K0Ut=2 — CUKPUTfc  START  LINE 

52  FORMAT ( Ic ) 

GO  TO  l 53, 54 , 84 , 500) ,KUDE 
500  K0DE=3 
oO  TO  53 

B4  CALL  STL2(XBl»Ybl,M) 

AT=3. 1415927*  !YU2(  1  )  +  Y  8  l  (  1)  )*ScRl  {  (Xb2(  l )  -  X  B  1  (  1)  )  *  *  2  + ( YB 2 (  1  )  - YB  1  (  I 

1  )  )  *  *  2 ) 

GO  1U  2 

53  REA 0(5, 1005)  {  X  P  (  J  )  »  Y  P  (  J  ) ,TH(J) , X  M  P ( J ) ,TPIJ)  ,PPlJ)  »  J  =  1  *  N  ) 

READ ( 5 , I0C5 } AT 

10C5  FORMAT ( 6c 13.0) 

DO  99  J  =  1 ,  N 

mK  1 TE ( 6, 79) XP { J) , YP{ J) , XMP ( J ) , J  H ( J ) , T  P ( J ) »  P  P ( J ) 

/9  (-URMAT (  IhO , 6 ( 3X,E15.8) ) 

99  TH  (  J  )  =  1  H  (  J  )  •  .  0  1  745329 
GO  10  2 
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LIST  OF  FORTRAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  BY  USING  THE  METHOD  OF  CHARACTERISTICS 


54  FE=FE*. 01745329 

RGi=(RT+YS)*COS<F£) 

RG2={RT-YS)*C0S(FE) 

Xl=(RT+YS)*SlN(Fh) 

X2=(RT-YS) *SIN(FE) 

AT  =  3. 14  1592  7* l RGI  +  RG2) *SORT( { X 1 -X2 ) *  *  2+ ( KG  1 -RG2 ) • *  2 ) 
CALL  SLRIN 
2  CP=GAM*R/GMl 
DO  60  J  =  1  *  N 
H=LP*TP( J)*G 
A=SURT (GMl*h) 

6C  VLP(J)=XMP(J)*A 
K'N 

DO  61  J=i,N 
FRX(  J)=XP(I<) 

FRYl J)=YP(KJ 
FRVt J)=VLP(K) 

FRT ( J ) =TP ( K ) 

FRP( J)aPP(K) 

FRTh{ J)=TH(K) 
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list  of  fortran  program 


plug  NO^ZLfc  ANALYSIS  BY  USING  THE  METHOD  OF  CHARACTERISTICS 

61  K=K-i 

F  L  X  (  n=XP(  l) 

FLYl  1 )  =  YP  (  I) 

F  L  P ( 1 ) =PP (  1) 

FUXI  1)=XP(M 
FUY ( 1)=YP(N) 

F  UP  I  I) =PP(N) 

GO  TU  [222 , 222,1 4)  , KODt 
222  M  =  N  +  N  -  1 
L  =  N  ♦  1 
J  =  0 

CALL  F  LDR  T  N  (l) 

M=  I 
L  =  N+l 

CALL  BPKTMI  ,  XB1 # YB1 , N l > 

FLX(NF)=XP{ I) 
r  L  Y l NF I =  YP  (  I) 

FLP(NF)=PPU) 

M  =  N-  i 
L~  2 
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LIST  OF  FORTRAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  BY  USING  THE  METHOD  OF  CHARACTERISTICS 


J*N 

CALL  FLDRTN  (  1 J 
M  =  N 

L=N+N-1 

CALL  BPRTN(2,XB2,Y62,N2) 
GO  TO  75 

74  CALL  BPRTM3»XB2»YB2»N2) 
NR  =  N+  1 
86  M=NR-l 
L  =  2 

J  =  NR- I 

CALL  FLDRTM3) 

M=  i 
L-2 

CALL  BPRTMi,XBI,Y81,Nl) 
FLX  J  NF ) =XP { 1) 

FL Y ( NF ) = YP (  l) 

F  L  P  i  NF  )  =  P  P  {  I) 

NF  =  NFf  I 
NR=NR+  l 
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LIST  OF  FORTKAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  BY  USING  THE  METHOD  OF  CHARACTERISTICS 


IFINR-IN  +  NU)  )  88 , 88 , 6 
88  I F ( J-666 186,6,86 
7S  FUX ( nF ) =XP ( N) 

FUY (NF J  =YP (N) 

FUP I NF ) =PP ( N ) 

NF=NF+  1 

IF ( NU ) 222 , 222  »  20 

1001  FORMAT (5E15.8/4E15. 8*312) 

1003  FORMAT ( 2t  1  5 . 8 ) 

1002  FORMAT ( 1 HO , / / / , 1 HO , 1  OX , 1H X , 1 7 X , 1 H Y , 1 7 X , 1 HM ,  13X  ,bH THETA,  13X, 
11HT, 17X, 1HP, 10X,3HI1R) 

20  KZ  ®  1 
NG=NF- 1 
NBP=N-2*NU 
NBZ=NBP/2 

I  F  (  (  N/2M2-N) 14,13,13 
1A  1F(2»NB2-NBP)22, 1,1 
13  IF  (2*NB2-NBP) 1 ,22,22 
l  NZ= ( NBP  +N- l ) / 2 
NY®  l 
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LIST  OF  FORTRAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  BY  USING  THE  METHOD  OF  CHARACTERISTICS 


GO  TO  3 

22  NZ= ( NBP  +  N ) 12 
NY  =  2 

3  M=N+N-1 
JL  =  N 
NB=  l 
L  =  N+1 
J  =  0 

CALL  FLDK  TN ( 2  I 
I  J  =  1 
1=0 
M*  l 

L  =  N  +  I  J 

CALL  BPRTNl 1 » X  B 1 »  Y  6  L  » N 1 ) 
A  FLX(NF )*XP ( I) 

FLY(NF)=YP( 1) 

FLP i NF ) =PP ( 1 ) 

NF=NF+ I 

A A  1 F ( J-666 ) A  5 »  6  »  A5 
A  5  J  =  L-1 
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LIST  OF  FORTRAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  8Y  USING  THE  METHOD  OF  CHARACTERISTICS 


L  =  2 
M*N  I 

I  F ( M-L ) 6  »  5 , b 
5  GO  TO  19,12) ,KZ 
9  IF(M-NZ)7,8,7 
8  KZ-2 

GO  TO  (11 ,12) ,NY 
II  CALL  FLDRIM2) 

J  =  0 

LXN+  l  J 
M  =  2*N  +  2*  I 
CALL  F LDRT N ( 2  ) 

NY  =  2 
M=  1 

L  =  N+  I  J 

CALL  BPRTNt 1,XBI,YB1,N1) 
GO  TO  A 

12  CALL  F  LDR  T  N ( 2 ) 

J  =  0 

L  =  N+ I J 
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LIST  OF  FORTRAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  BY  USING  THE  METHOD  OF  CHARACT t;R I  ST  I CS 


M=2*N+2*I-1 
CALL  FLDRTM2) 

M=  1 

LCEN+  I  J 

CALL  BPRTN(1,XBI,YB1#NI) 

I  J  =  1  J—  1 
1  =  1-1 
GO  TO  A 

7  CALL  F  LDR T  N ( 2 ) 

J  =  0 

L  =  N+  I  J 
M  =  2  *N  +  2*  I 
CALL  FLDRTM2) 

M=  1 

L  =  N+  l  J 

CALL  BPRTNt 1 ,  XB L , YB 1 , N 1 ) 
I J  =  1  J  +  1 
1  =  1  +  1 
GO  TO  A 
6  SUMP  =  0 . 
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LIST  OF  FORTRAN  PROGRAM 


plug  nozzle 


ANALYSIS  BY  USING  THE  METHOD  OF  CHARACTERISTICS 


writer, 3051 

305  FORMAT (1H0/1H0»3AHTHRUST  DISTRIBUTION  ALONG  THE  PLUG) 

IFIKODE  .EQ.  3)  GO  TO  207 
I  =NG 
NG=NG- 1 
DO  64  K= I »  NG 
1  =  1-1 

P 1 2= ( FUP ( I ) +FUP { I+L) )/2. 

Al2.J.UlS«7.|FUVimFUV«I*l)).S8RIUFUXCI>-FU*n*l))..2*tFUVU. 

L-FUY( l  +  l)  ) **2  5 

FEI?  =  A1  ANMFUXI  I  J-FUXI  I  +  l)  )  /  (  FUY  (  I  ♦  i  J-FUY  (  I  )  )  ) 
PI=P12*A12*C0S(FE12) 

IFIFUYII+L)"FUY(I))71»76»70 
76  Pl=0. 

GO  TO  70 
71  P I =-P  I 
70  SUMP  =  SUMP  +  P l 
64  CONTINUE 
207  SUMM  =  0  * 

$UMV=0. 
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list  of  fortran  program 


plug  NOZZLE  ANALYSIS  BY  USING  THE  METHOD  OF  CHARACTERISTICS 


NZ=N-1 

DO  91  1  =  1  *  N/ 

PL2MFRPI  l  )+FRP(  1  +  1)  )/2. 

T 1 2= { FRT ( I  )  +FR T (  1  +  1)  )/2. 

R012=P12/(R*T12) 

4l2  =  3.1415927MH<m)*FKY(lU))«SQRII(KRXII)-FRX(m))-2MFRY(U 

1-FRYl I+L) )+*2) 

TH12=(FRTH(  I  )+FRTH(  I  ■*  l)  )/2. 

FE12-ATANI  ( FRX {  I  ) -FKX (  I  +  1 )  ) / ( FRY ( 1 +1 )-FRY(I) ) ) 

V  l  2  = ( FRV (  1  )  +PRV (  1  +  1)  )/2. 

SQ=FEl2-TH12 
VM=R012*VL2*A12*C0S (SQ) 

VMQM=VM/G*V12*C0SI TH12 ) 

VMQMP3Pl2*Al2«Ct)S(Fti2) 

VM(J  =  VMUM  +  VMUMP 
SUMM=SUMM+VMO 
SUMV=SUMV+VM 
91  CONTINUE 

SUKM=SUMP+SUMM 
CF I =  S U M M / (  P*AT  ) 
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LIST  OF  FORTRAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  BY  USING  THE  METHOD  OF  CHARACTERISTICS 


WRlTfc(6,96)SUMM,CF l 
WRITE ( 6, 400 ) SUMV 

400  FORMAT {  IHO, I5HMASS  FLOW  RATE=£15.8) 

9  6  FORMAT  (  IHO  ,  !>HSUMM=t  15. 8 , 3X  ,  4HCF  1  =  E  1 5 . 8  ) 

SUMP=0. 

NB=NF- 2 
00  92  1=1, NB 
IF(SUMP)602»602,60L 
602  CONTINUE 

IF (FLX  l  1)-FLX(  14  l)  )600,92,92 

600  IF ( I .EU. L )CU  TO  6G  L 
P12=(FLP(1)+FLP( 141)1/2. 

AL2  =  3.  L415927*  (FLY  ID+FLYt  I  +  I)  )  »SQRT  (  (FLXm-FLX{  1  +  1)  )  »*2 
l  +  IFLYI  1  )  —  FLY (  l4i) )**2) 

FE  12  =  ATANUFLXin-FLXI  1*1)  )  /  t  FLY  11*1 ) -FLY  (  I  )  )) 
PI=P12*AL2*C0S(FEI2) 

1FIFLYI 1)-FLY( I+L) >72,73,73 

601  P12=lFLP{ I )+FLP{ 1*1) )/2. 

A  1 2  =  3  .  L4L5927MFLYI  I)+FLY(  I  + 1 )  )  •  SORT  (  I  F  LX  C  I  )  — F  LX  {  I  4- 1  }  )  *«  2 
14 ( FLY(  I  )-FLY(  I 4J  ) ) »*2 ) 
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FE  12  =  AT  AN  I  (I- IX  (  l  )-FLXt  I  *1 )  )  /  I  FLY  l  l  +  l  J-FLYl  I  ))  ) 

PI=P12*A12*C0S(F612) 

IF(FLY(I)-FLY(I+L))72,73.73 

72  P  l  =-P I 

73  SUMP=SUMP*PI 
TOT=SUMM+SUMP 
CF=TUr/(P»AT) 

WRITEI6»94)FLX(I  +  IJ  *FLYII+l)»TUI  »CF 

92  CONTINUE 

9  A  FORMAT  (  1H0»?HX  =  £1. 5«8t3X»2HY  =  Fl*j.8»3X»2HT  =  EI‘j.8f3X»3HCF  =  L:  IS. 8) 
AX=AT»144. 

WRlTE(b,93)AX 

93  FORMAT ( 1H0,3HA1 =ELB. 8) 

DUMMY  =  XMP (  l ) 

CALL  CGNVR  ( I » DUMMY » W  l ) 

THA-THl 1 ) +WI 
IFI THA)20l,202,202 
201  WRITE16,200) 

PB-PA 
GO  TO  203 
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200  FORMAT! IHO, 10X,34HFLOW  BREAK-AwAY  FROM  UPSTREAM  W ALL/ IHO , 10X , 9HSE T 
I  P  B  =  P A ) 

202  CALL  CUNVK (2 , XMRED, THA ) 

PRtQ=L./(  ( l.+GMl/2.»XMKED**2) *  * (GAM/GM1  ) ) 

P  0  A  =  { 1.+GMI/2. *  X  M  P  ! 1 ) *  »  2 ) »  * (GAM/GM1) 

CALL  B PR S 

PB  =  T  ABLE  1  ( PBP1 ,XM1 ,XMREO,NC) 

P2PR=PB 

P2P l=P2PR*PUA*PRED 
I  F  (  P  ?P  I-  1 .  )?04,205,205 

205  WR 1 T t ( 6 , 2 06 ) 

PB  =  PA 

GO  TU  203 

206  FORMAT ( LHO, i 0 X , 2 2H THEORY  BECOMES  I N V AL I  0/ 1H0 , l OX , 9H SE T  PB  =  PA) 

20 A  PB  =  PP( 1 ) *P2P  I 

203  CONTINUE 

TVAC=T0T+PB*3.1AI5927*FLY(NF-1)*»2 
CFVAC=TVAC/IP*A1 } 

AE=3.1A15927»YB2(N2)“*2 

THRUST=TVAC-PA*AE 
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CFREAL *THRUST / { P*AT) 

PB*PB/144. 

WRITEC6,6lll)PB,TVAC,CFVAC,THRUST,CFREAL 
6111  FORMAT(1HO,3HPB=E15.8/1HO,5HTVAC*E15.8/1HO,6HCFVAC=E15.8/1HO,7HTHR 
lUST*E15.8/lH0f7HCFREAL“E15.8) 

SI  kR I TE ( 6  *  2000 ) 

2000  F0RMAT(1H0*///,20X,10HEND  OF  JOB) 

GO  TO  78 

END 

SUBROUTINE  SLRTN 

DIMENSION  YP(400) .XPI400) ,TH(400) ,XMP(400> ,TP(400) ,PP(400) , 
1RXM{200)*RTH{200)»RTP(200),RPP(200) ,VLP(400) 

COMMON  YP,XP,TH*XMP,TP,PP,RXM,RTH,RTP, RPP , VLP, ROS , YS, GAM, GM1, G, 
iXM,N,P,T,R,L,M,J,N2,XX82,YYB2,NU»KNT,GPl,FE,RT 
1«  PA 

IFIN-*  10)50*51*51 

50  XN*N 

GO  TO  52 

51  XN»N-4 


52  CONTINUE 
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PM^P*(2./GPI)»*(GAM/GMI) 

TS  T  =  T* ( 2 . /GP 1 ) 

A  =  SORT  {  I./  (GP1*RGS»Y$)  ) 

EPP-YS/ 6. *SQRT (GPI*YS/RUS) 

FEU-ATAN [ tPP/KT  ) 

F I=FE+FEO 

RTG  =  SQRT (LPP*tPP*Kl*KT  ) 

HH  =  RTO*SIN( F  I  ) 

HK=RTO*COS {FI) 

XMt  S  =  SUK1  (  (GPl/2.»XM«XM)/{l.+GMl/2.*XM*xM)  ) 
XPP=(XMf.S-  l  .  )  /A 
PHA-ARSIN! (LPP+XPP)/ROS) 

YL  =  YS  +  (R(jS-R(jS*CQS  (PHA)  ) 

D  Y  L  -  2 « »  YL  / lXN-1. ) 12. 

DO  l  I  =  1 »  N 
IF  (  N-  1  0  )  A  A  ,  9  »  9 
9  IF {  1-6) A , 2 ,3 
3  IF (  I  —  ( N  ~  3  )  )6,5,  7 
2  UYL=2.*DYL 
GO  TO  6 


CHARACTERISTICS 
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5  DYL=DYL/2. 

GO  TO  7 

6  YPP=FLOAT l 1-3) *DYL-YL 
GO  TO  8 

7  YPP=FLOAT( I+N-10) »0YL-YL 
GO  TO  8 

44  YPP=FLOAT{ I-I ) «DYL»2.-YL 
GO  TO  8 

4  YPP=FLOAT (1-1 ) »DYL-YL 

8  U=A*XPP+GP 1/2 • *A*A*YPP«YPP 
V=A*A*GPI*(XPP*YPP+GP1/6.*A*YPP»*3  ) 

XM  S  =  SORT l { I.+U)*«2  +  V*V) 

THX=ATAN( V/( I.+U) ) 

XMP (  I)=S0RT(2./GPl*XMS*XMS/l l.-GMl/GPl*XMS*XMS) ) 
■PP(I)*PST/{ {2./GPl)**(GAM/GMl)*( l.+GMl/2.»XMP( I )**2) 
I** (GAM/GM1 ) ) 

TP(I)  —  TST/( l I . +GM1 /2. • XMP { I ) **2 ) * { 2. /GP1 ) ) 

XPI  I }  =  XPP*COS (-F I )  - YP  P  *  S I N ( -F I ) +  HH 
YP (  I  )=YPP*COS{-FI )+XPP*SlN{-Fl  )+HK 
THU  )  =  THX-F  i 
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THH=TH( I )*57. 29578 
QX=XP { I ) *  1 2 . 

QY=YP( I ) *  1 2 . 

UP=PP ( I ) / 1AA. 

I  WRITE (6,  L02  )  QX  ,  QY  *  XMP  l  I )  ,THH,TP(  l )  ,CP 
102  FORMAT { 1H0,6( 3X.E15.8)) 

RETURN 

ENO 

SUBROUTINE  S T L2 ( X B , Y B , N 1 ) 

DIMENSION  YP(AOO) , XP ( 400 ) , TH ( 400 ) , XMP (400) , TP ( AGO ) ♦ PP t AOO ) , 
1RXMI200) , RTH ( 200 ),RTP 1200) ,RPP(200) ,VLP(400) 

COMMON  YPfXP»TH»XMP»TP»PP»RXMfRTH»RTP»RPP»VLP»ROSfYSfGAM»GMl,G» 
IXM,N,P»T.R,L*M,  J,N2»XXB2»YYIi?,NU,KNT,GPl,FE,RT 

1 »  P  A 

UIMtNSION  X6(  100) , YB( 100) 

P I =3  «  1 A  1 592  7 
READ15,  IDKODt 
GO  TO  112,13) , KODE 

II  FORMAT ( 12) 

13  THl  l)  =  FE 
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FE=FE-ARSIN( l./XM) 

GO  TO  16 

12  THST=FE+PI/2. 

FE=FE+(Pl/2.-AKSIN( l./XM) + SORT (GP1/UM1 ) *ATAN ( SURF { GM1/GP 1  *  t  XM»XM-  1 
l . )  ) )-AT AN ( SORT ( XM*XM-l.  )  )  ) 

16  JJ=1 

TSC=SIN(P1+FE)/C0S(P!+FE} 

1  YD=YB I  JJ+1  )  -  YB  {  J  J ) 

XD=XB(  JJM  )  —  X b  (  JJ) 

XA=  1 .  /  {  YD/  XD-  TSCJM  YYB2- YB { J J ) *YO/XD*XB (JJ)  -XXB2*TSC  ) 

IFIXA-XBI JJ)  ) 3  »  2  t  2 

2  IF(XA-XtH  JJ+1)  )4,4,‘j 

5  JJ=JJ+1 
GO  TO  1 

3  WRITE{6,6) XA,XB( JJ) ,XBI JJ+1  ) 

6  FORMAT l 1H0,3E15. 8) 

STUP 

A  YP  l  1  )=YB ( J J)  +  ( XA-XBt JJ) ) /XD*YD 
XPl  1 ) =  X  A 


GC  TO  114*10) » KUDE 
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1 A  TH  t 1)  =  THST  +  SCR I (GPI/GMi ) *ATAN(SURT (GMI/GPi* (XM*#2-1 . ) ) )-ATAN(  SORT  ( 
1XM**2-1. ) ) 

15  XMP ( 1)=XM 

PP(l)=P/( ( l.+GMl/2.*XM»*2) *  * ( GAM/ GM I )  ) 

TP(l)=T/(l.+GMl/2.*XM**2) 

THP=TH ( 1) *57.29578 
QX  =  XP ( 1 ) *  1 2 . 

OY  =  YP( 1 ) *  1 2 . 

GP  =  PP  (  D/144. 

WRIT£(6* IO)0XtUY«XMP(l) »  THP  »  TP ( I) ,QP 
XN  =  N 

XN  =  XN-  1 . 

0X=(XXB2-XA)/XN 
DO  9  MM  =  2  »  N 
XP{MM)=XP(MM~1)+DX 

YP(MM}=YP(MM-1)  +  (YYB2-YP(  I))/(XXB2-XPU))*l)X 
XMP ( MM ) =XM 
TH ( MM) =TH {  I) 

TP (MM) =TP (  l ) 

PP ( MM) =PP  t  1  ) 
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QX=XP(MM) *12. 

QY  =  YP l MM ) *  12 . 

9  WRITE(6»i0)QX»UY»XMfTHP,TP{ l) ,  UP 
10  FORMAT l lH0t6(3X.E15.8) ) 

RETURN 

END 

SUBROUTINE  &PRTN(KUD£,XfY,NX) 

DIMENSION  YPJ400) , XP1400) , TH(4uO) ,XMP{400) , T P ( 4 00 ) , PP l 4 00 ) , 
1RXM1200) ,RTH(200) ,RTP(200) ,RPP(200) ,VLP(400) 

DIMENSION  XI 100) , Y (  IOC) 

COMMON  YP, XP,TH,XMPrTP,PP,RXMfKTH,RlP,RPP,VLP,ROS»YS,GAM,GMl,G, 
lXM,N.P»T»R»LtM, J»N2»XXU2tYYB2.NU,KNT,GPltFE,RT 
1  *  PA 

IF (KODE  .EQ.  3)  GO  TO  65 
I  TR=  1 

WR1TE16, 1004) 

1004  FORMAT ( 1H0.40X, 18HB0DY  POINT  ROUTINE) 

CP^GAM»R/GM1 
H1~-CP*TP!L)*G 
A 1  =  SQR  T ( GM 1 »H 1 ) 
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V1=XMP(L) *A1 
B1=ARSIN(A1/V1) 

H=Hl+Vl*Vl/2. 

P$=P*(2./GP1)**(GAM/GM1) 

TS=T*(2./GPL) 

Sl=CP*ALOG( (TPIL)/TS)/{{PP(L)/PS  )**IGMI/GAM) ) ) 
SB=CP*ALOG((TP(M)/TS)/UPP(M)/PS  )  **  t  GMl/GAM 

DO  1  1 -  1 »  N  X 
K-  I 

IF  IX( I )-XP (L)  )  If I»2 

1  CONTINUE 

20  GO  TO  ( 60t 61  *  52) »KUOE 

60  KE Y=  1 

GO  TO  23 

61  KE  Y  =  2 

GO  TO  23 
23  K=NX 

TH3=ATANl (Y(K)-Y(K-l) J/CX(K)-X(K-l) ) ) 

B3= ARS IN l L./XMPIM) ) 

A3=SURT (GMl*CP*TP(M)*G) 


CHARACTERISTICS 
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\ 

V3=XMP l M ) * A3 
T3=TP l M ) 

TH4^TH(M1 
58  CONTINUE 
3  7  GO  TO  I  34 ♦ 35 ) ,KtY 
35  C 1  =  ( TH3  +  B3*  TH4  +  B4I/2. 

GO  TO  43 

34  CL=(TH3-B3+TH4-B4)/2. 

43  Ci=SIN(CIJ /COS(Ci) 

8=1  {Y(K)  —  YP(L)  )  —  (X (K)  —  XPCL)  ) * C i ) / t (YP(M)-YP(L)  )— (XP(M)-XP(L)  )*CL) 
X4=XP(L)+B*{XP(M)-XP(L) ) 

Y4=YPl L) +B»( YP(K)-YPI L) ) 

XM4=XMP (L ) +  B* ( XMP ( M)-XMP( L ) ) 

TH4=TH(L)+B*( THIM)-TH(L) ) 

T4=TP(L)+B*(TP(M)-TP(L) ) 

P4=PPIL)+B»(PPCM)-PPCL) ) 

S4  =  S  I  +  B* l SB-S 1 ) 

TH44=TH4*57.2S578 
A4=SQRT (GMl*CP*T4»G) 

V4=XM4*A4 
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B A  =  ARS I N (  l./XMA) 

GO  TO  { AS, 46) ♦ KE Y 
A 6  C3=COS(  {TH3  +  B3+-IHA  +  BA)/2.  ) 

GO  TO  A  7 

A 5  C3  =  COS ( l TH3-B3  +  THA-BA) /2.  ) 

A  7  C  l  = { B3  +  BA ) /2 . 

C2  =  COS ( C 1 ) 

C  1  =  S  I N  (  C  i  ) 

GO  TO  ( A 8  »  A 9 )  ,  KEY 

A  9  V3-VA+ (  (V3  +  VA)/2.*C1/C2)*( (TH3-IHA  +  C1*SIN( (TH3  +  THAJ/2. )/(  (Y(K) 
1  +  YA) /2. »C3 ) » ( X(K)-aA)-( ( T3  +  TA ) /2. )/( ( ( A3  +  AA) /2.  ) «*2) »CI*C2  * 

1 ( SU-SA ) *G ) ) 

GO  TO  50 

A  8  V3  =  VA+( (V3  +  VA)/2.*CL/C2)*( (THA-TH3  +  Cl*SIN( ( IH3  +  THAJ/2. )/( (Y(K) 
l+YA)/2.*C3)*(X!K)-XA)-{(T3+TA)/^.)/{ ( (A3+AA)/2.)**2)*C1*C2  * 

1  f  SB-SA ) «G ) ) 

50  H 3  =  H~ . 5  *  V  3  *  V  3 
A3  =  SURT (GMI*H3  ) 

B3=ARS1N(A3/V3) 

XMP3=V3/A3 
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T3=H3/ ( CP*G ) 

P3=PPIM)*(T3/TP(M) )»*{GAM/GMI) 

THP=TH3»57. 29578 
IF ( ABS (0-BP )-. OOOOO 1)56,56,57 
57  6P=8 

GO  TO  58 
56  WRITE(6,205) 

205  FORMAT ( IHO, 12HINSEKT  POINT) 

UX=X4*i2. 

U Y- Y4* 12 . 

QP=P4/144. 

WRIT£{6,I006)t!X,QYtXM4,TH44,T4,CP 

WRITEI6.2C6) 

206  FORMAT ( IHO, I2HC0RNtR  POINT) 

UX=X(K)*I2. 

UY-Y ( K ) *  12 . 

QP=P3/144. 

WRITE! 6, 1006)QX,QY,XMP3,THP,T3,CP 

VLP(M)=V3 

XP ( M ) *  X ( K ) 


58 


LIST  OF  FORTRAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  BY  USING  THE  METHOD  OF  CHARACTERISTICS 

YP { M ) =  Y ( K ) 

XMP { M ) =XMP  3 
TH(M)=TH3 
T  P  (  M  )  =  T  3 
P  P  (  H  )  *  P  3 

GO  TO  (51,52) ,K£Y 

51  J  =  666 
RETURN 

65  TH3=TH(N) 

XMP3=XMP(N) 

NZ  =  N 
P  3  =  PP ( N ) 

T  3  =  T  P  l  N ) 

K  =  NX 

X ( K ) =XP ( N ) 

Y ( K ) =YP ( N ) 

52  READ(5,1007)NU,PA,KKD 
IF (KKD.EU.OJGO  TO  401 
PA=PA* 144 . 

401  WRITE(6,207) 
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207  FORMAT  (  JLHO  ,  29HR I GHT  RUNNING  CHARACTERISTICS) 

1007  FORMAT ( I2.E15.8, 12) 

C3= ( GAM+ 1 • )/GMl 
CW./C3 

TERM=TH3-SGRT  C  C  3 ) * AT AN { SQR T{ C4 *  I XMP 3  *  *  2- 1 . ) ) ) ♦ AT  AN ( SORT  { XMP 3* *2-1. 

in 

XKt  =  SQRT(2./GMl»{ (  l.+GMl/2.*XMP3**2)/( ( PA/P3 ) *  * ( GM 1 /GAM ) )-l. ) 

1) 

XNU=NU 

DM=(XME-XMP3)/(XNU-1.  ) 

XM  =  XMP  3 

53  DO  54  II*1»NU 
RXM ( I  I ) =XM 

RTH( II }’TERM+SQRTtC3)*ATAN{SQRHC4*{XM*XM-l. ) ) )  - AT AN ( SORT ( XM*XM- 1 . 
L)  ) 

RTP( II )=T3/(1.+GMI/2.*XM*XM)*C1.+GM1/2.*XMP3*»2) 

RPPI  I  I )~P3»  C ( l.+GMl/2.*XMP3**2) / ( 1.+GM1/2. *XM*XM) ) 
l»» (GAM/GM1  ) 

RTHP=RTH ( I  I ) *57. 29  578 
XM=XM+DM 
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YB  I  =  Y { J ) 

YB2=Y(K) 

XB=XB 1 
BB  =  B  I 
VB  =  V  1 
TB*  TP l L  ) 

AB  =  A  1 
KKNT=0 

II  KKNT  =  KKNT ♦  1 
KCNT  =  0 

IF(KKNT-50)lll,lllf 133 
133  WRJTE(6,134)XBP,XB,XB1,XB2 
.13 A  FORMAT  (1H0,4E15.8) 

GO  TO  13 

ill  TH{Mj=ATAN{ (YB2-YBI)/(X82-XB1  )  ) 
22  GO  TO  (33, 44) , KODE 
33  Zl=(TH(L)-Bl+TH(M)-bB) /2. 

Z9=TH< L)-TH(M3 
GO  TO  55 

44  Zl=( TH(L)+Bl+TH(M)+bB)/2. 


62 


LIST  OF  FORTRAN  PROGRAM 


PLUG  NOZZLE  ANALYSIS  BY  USING  THE  METHOD  OF  CHARACTERISTICS 


Z9*TH( MJ-TH { L ) 

GO  TO  55 
55  Z2*C0S(ZL) 

KCNT  =  KCNT  + 1 

Z10  =  SIN{ (Bl+BBJ/2.  ) 

Z3=SIN(Z1)/Z2 
Z4=SIN(THIM) }/COS(TH(M) ) 

Z5*S  I N  ( (TH(L)+TH(M) )/2.  ) 

Z7*(Bl+BB)/2. 

Z8*C0S ( Z 7  ) 

Z6=»SIN  ( Z7 ) 

XBP=(YP{  U-~YB1-XP<  LJ»Z3  +  XB1*Z4)/{Z4-Z3) 

YP(M)=Y81+(XBP— XBl)»Z4 

VBP  =  VI+( IVl+VB)/2.*Z6/Z8)*(Z9+(Z5»Z10/( (YPILJ+YPIM) )/2.*Z2) ) ) *  ( 
IXBP-XP ( L ) )-((TP(L)+TBJ/2./{( { A I +AB ) /2 . ) » *2 ) *Z6*Z8 ) * { SB-S I ) *G 
12  HBP  =  H— VBP*  VBP/2. 

AB  =SQRT { GM1*HBP ) 

XMP (M)=VBP/AB 
T8  =GM1*HBP/(GAM»R*G) 

IF < XMP 1997,998,998 
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997  WRITE16,1006)XMPIM) 

STOP 

998  CONTINUE 

-IF  ( ABS  (  (XBP-XB)/XB)-.000001)3.3,4 
4  BB*ARSINU./XMP(M)  ) 

XB-XBP 

ITR-ITR+I 

VB-VBP 

IF.(KCNT-50)22,22,333 
333  WRITEI6, L34 ) XBP , XB , X8 l , XB2 
3  IF(XB1-XBP)6,13,5 
6  IFIXBP-XB2) 13, 13,9 
9  XB laXB2 
YB I* YB2 
J*J  +  1 
K  =  K+ 1 

I Fj{ K-NX)  21 ,21 ,200 
200  TH ( M ) “XTH3 
YP { M ) *XY3 
XMP ( M ) *XXMX3 
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GO  TO  20 
21  XB2=X( K } 

YB2*Y4K) 

GO  TO  11 
5  XB2=XB 1 
Y62*YB1 
J*=  J-  L 
K»K-1 
X8 1=X { J ) 

YBI=Y{ J) 

IF ( J ) 20 1  20  *  1 1 

13  THB2=TH(M)*57. 29578 

PPIM)=PP{M)*(TB/TP(M) ) *« ( GAM/GMI ) 

TP (M)«TB 
QX  =  XBP  *  12 • 

QY=»YP(M)*12. 

QP=PP(M)/1AA. 

WRITE! 6, 1006)QX,QY,XMP(M) ,THB2,TPIM) ,QP, ITR 
VL  P ( M ) = VBP 
XP (M)*XBP 
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RETURN 

1006  FORMAT(1HO,6(3X,E15.8) ,15) 

END 

SUBROUTINE  FLDRTN  (IZ) 

DIMENSION  YP ( 400 ) ,XP(400) , TH ( 400 ) , XMP { 400 ) , TP ( 400 ) , PP ( 400  )  , 

1RXM ( 200 ) , RTH ( 200 ) ,RTP{200) , RP P { 200 ) , VLP ( 400 ) 

DIMENSION  H(3) »A(3) , V ( 3 ) , B { 3 ) ,S(3) 

COMMON  YP,XP,TH,XMP,TP,PP,RXM,RTH,RTP|RPP,VLP»R0S,YS,GAM,GM1,G, 
lXM,N*P,T,R,L,Mt J,N2,XXB2,YYB2,NU,KNT,GPI,FE,RT 
1 ,  PA 

WR ITE ( 6 , 2 ) 

2  F0RMAT(1H0,40X,13HFIELD  ROUTINE) 

MS*1 

CP*GAM*R/GM1 

GO  TO  (23,23,24) , IZ 

23  1 1 *L 

GO  TO  25 

24  I  I *M 

25  DO  10  I  J*L  » M 
I  TR=  i 
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GO  TO  (32,32,18) ,IZ 
18  J=J-1 
GO  TO  11 

32  J*J+i 

GO  TO  33 

33  GO  TO  ( 11, 12)  ,  IZ 

12  IF(II-M) LI, 13, 13 

13  IF(KNT-NU) 14, 14, 11 

14  M5*2 
ST=TP( J+l) 

SP=PP( J+l) 

SM  =  XMP ( J  + 1 ) 

SH=TH( J+l) 

SX=XP( J+l) 

SY=YP( J+l) 

TP( J+1)=RTP(KNT) 

PP ( J+l )=RPP(KNT) 
XMP(J+1)=RXM(KNT) 
TH( J+l )=RTH(KNT) 

XP( J+l )=XXB2 
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YP ( J+ 1 ) =YYB2 
11  DO  8  1*1 i 2 

GO  TO  (19, 19, 20), IZ 

19  J1*J+I-1 
GO  TO  21 

20  J.l*J  +  2*  1-2 

21  H ( I)*CP*TP( J1)*G 

A ( I )=SQRT(GM1*H( I ) ) 

V( I)*XMP( J 1 } * A  C I ) 

B(I)*ARSIN(A( I)/V( I)) 

PT*P*(2./GP1)»MGAM/GM1) 

TT*T*(2./GP1) 

S(  I)*CP*ALOG( (TP ( Jl)/TT)/( ( PP l J 1 ) /PT ) ** ( GM 1 /GAM ) ) ) 
8  CONTINUE 

W*H(l)+Vm*Vll>/2. 

TP  C 1 1 ) =TP ( J) 

B  (  3 )  *  B  (  1 ) 

TH( II) =TH  C J) 

S ( 3 ) *S ( 1 ) 

A ( 3 ) *A ( 1 ) 
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V I  3 ) =V ( 1 ) 

4  GO  TO  (26,26,27) ,  IZ 

26  JP*J+1 
GO  TO  28 

27  JP*J+2 

28  Z1S(TH(J)+B( 1 ) +TH { 1 1 ) +B ( 3 ) )/2. 
Z2=(TH{JP)-B(2)+THlII)-Bl3))/2. 
24= ( B ( 1 ) +B ( 3 ) ) /2. 

Z5=(B(2)+B(3) )/2. 

Z6=( V( 1)+V(3) ) /2. 

Z7= ( V ( 2 ) +V { 3) )/2. 

Z I2=COS  l  Z 1 ) 

ZI3=COS(Z2) 

ZI6=C0S(Z4) 

Z17=COS(Z5) 

5  FORMAT ( 1H0,6( 3X,E15.8) , 15) 

Z8  =  S  IN(Z1)/Z12 
Z9=SIN(Z2) /ZI3 

Z 1 Q=S I N ( Z4 ) 


ZI1=SIN(Z5) 
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ZIW16/Z10 
Z15*Z17/ZJL1 
Z18  = ( TP { J ) +TP { 1 1 ) )/2. 

219*(TH( J)+TH(II))/2. 

Z2Q=(TH(  JPJ+THUI)  )/2. 

Z2132. «Z  18 

Z22=*((A(l)+A(3))/2.  )**2 

XP( II)=(XP( J)+1./Z8*(YP( JP)-YP( J)-XP( JP)*Z9))/( 1.-Z9/Z8) 

Z25=XP  C I  I )-XP ( JP) 

Z26=XP ( II)-XP(J) 

YP(II)*YP( JP)+Z9*Z25 
Z23= ( YP { J ) +YP ( 1 1 ) ) /2. 

Z24= ( YP ( JP ) +YP (  I  I )  )/2. 

S(3)>S(l)+( ( S { 2 )-S  £  1 ) )*Z26*(Z10/ZI2) ) / ( Z26*Z 10/Z I2  +  Z2 5*Z I l/Z 13 ) 

VI 3)=I./(Z1A/Z6  +  Z15/Z7)»(TH( JP)-TH{ J ) +ZI4/Z6*V ( 1 1 +Z 15/Z7»V ( 2 )  + 
1Z10*SIN(Z19)/(Z23»Z12)*Z26+Z11*SIN(Z20)/(Z2W13  W25-Z18/IZ22)  *Z 
210*Z16*ISI3)-S(1)) *G- (Z18/Z22)*Zll*Z17*(S(3)-S(2))«G) 

TH3P=TH( J)+(V(3)-V(I) >  *  I Z14/Z6 ) - I Z 10*SINl Z 19 } ) / i Z2  3*Z 12 ) *Z2  6  +  Z 1 8/ 
1Z22*Z10*ZI6*(S(3)-S(1) J *G 
H(3)=W-VI3)*V(3)/2. 
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A( 3)=SQRT(GM1*H(3) ) 

TP{ II)=GMI*H(3)/CGAM*R*G) 

IF( I TR-50) 67,67,68 
68  WRITE(6,5)TH3P,TH{ Ili 
GO  TO  6 

67  I F ( I TR- 1)7,7,66 

66  IF(ABS(TH3P-TH( II) OOOOO L ) 6 , 6 , 7 
7  B(3)=ARSIN(A{3)/V(3) ) 

TH{ 1 1 J  =  T H 3 P 
I TR= I TR+ 1 
GO  TO  4 
6  TH { 1 1) =TH3P 
VLPt  m=V(3) 

THPP=*TH ( I  I ) *5T . 29578 

PP( 1 1 J  =  PP( J)*t  CH(3)/H(1) J**(GAM/GM1) ) / ( EXP ( l S t 3)-S ( I) ) /R) ) 
XMP { II)=V13)/A(3) 

QX*XP(II)*12. 

QY* YP ( I  I) *12. 

QP  =  PP { I  I ) /  144. 

WRITE{6t5)QX,QY,XMPlII ) ,THPP,TPC II ) ,QP*ITR 
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GO  TO  {30,30,31) ,IZ 

30  11*11*1 
GO  TO  10 

31  11=11-1 
10  CONTINUE 

GO  TO  (17,16) , MS 

16  TP { J  +  l ) =ST 
PP(JU)*SP 
XMP( J+1)=SM 
TH( J+l )=SH 
XP{ J+1)=SX 
YP(J+l)=SY 
KNT=KNT+i 

17  RETURN 
END 

SUBROUTINE  BPRS 

DIMENSION  YP(AOO) , XP { 400 ) , TH ( 400 ) , XMP ( 400 ) , TP ( 400 ) , PP ( 400) , 
1RXM(200) ,RTH(200),RTP(200) ,RPP1200) ,VLP(400) 

DIMENSION  ZC( 100), ZJ( 100) ,XM1( 100) ,PBP 1(100) 

COMMON  YP,XP,TH,XMP,TP,PP,RXM,RTH,RTP,RPP,VLP,ROS,YS,GAMsGMl,G, 
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1XM,N»P,T,R» L,M, J,N2,XXB2,YYB2,NU,KNT,GPI,FE,RT 
I, PA 

COMMON ZC,ZJ»XMlfPBPltNQ 
ON*.  I 

K«1 

XM2=I.5 

2  CZ=XM2»*2/(2./GMl+XM2»*2) 

QJ=TABLEl(ZJ,ZC,CZfNQ> 

CD=«Q  J*QJ*CZ 

zz*i.y 1 1 i.-cd) *  » igam/gmi ) ) 

TZZ=U5.*IZZ-1.)  )/(7.»XM2**2-5.*CZZ-I.  )  )  > »*2* I ( 7 . *XM2»*2- t 6 . »ZZ  + 1 . 
U )/(6.*ZZ+I.  )  ) 

T 34= AT AN ( SQRT(TZZ) ) 

WI»SQRT1GP1/GMI)»ATAN(SQRT(GML/GPI*(XM2«»2-1- ) ) )-ATAN( SQRTIXM2 
1*«2-1. ) )  -T  34 
CALL  CON VR l 2  »  XM »  W l ) 

P1P0=(  l.+GMl/2„*XM**2) 

P0P2= ( l.+GMl/2.*XM2**2) 

PBPl(K)=(PlPO/PQP2)**(GAM/GMl) 

XM1 CK)=XM 
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XM2=XM2+DM 
K*K  +  1 

IF ( K~100 ) 7  »  7  » I 
7  IFIXM-6. )2,1»I 
1  NQ*K-l 
RETURN 
END 

SUBROUTINE  CON VR { KODE , XM , ANGLE ) 

DIMENSION  YP(400) ,XP(400) ,TH(400) .XMPI400) , T P ( 400 ) , PP { 4 00 ) , 
1RXM(200) ,RTH(200) ,RTP(200) .RPP1200) ,VLP(400) 

COMMON  YP,XP,TH,XMP,TP,PP,RXM,RTH,RTP,RPPTVLP»ROS,YS,GAM,GMl,G, 
IXZ,N,P,T,R,L,M, J,N2,XXB2, YYB2 ,NU, KNT , GP1 ,FE,RT 
1 » P  A 

GO  TO  (1,2), KODE 

C  K0DE=1 — INPUT  M,  COMPUTE  ANGLE 

C  K0DE*2— INPUT  ANGLE .COMPUTE  M 

1  ANGLE* SORT ( GP1/GM1 )*ATAN( SORT { ( GMI /GP 1 ) * { XM*XM-1 . ) ) ) 

1-AT.*M  (  SORT  ( XM«XM-1.  ) ) 

RETURN 

2  XM* 10. 
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J-0 

K£Y=0 

DXM=*1. 

55  IF ( J-50 ) 5 , 13,13 

5  FMI»SQRT{GPl/GMi)*ATANlSQRT{GMI/GPI»(XM«XM-i. ) ) J-ATAN t SORT ( XM*XM 

1-1.  J  ) 

TEST=FM I-ANGLE 
I  F.(  KE  Y  )  4 , 4 , 3 
4  XM=»XM-DXM 

IF(TEST)8, 13,9 
9  KE Y= 1 
GO  TO  5 
8  KEY 3 2 
GO  TO  5 

3  GO  TO  (6, 7), KEY 

6  I F { T  ES  T  3 10,13,11 
11  XM=XM-DXM 

J  =  J+1 

IF ( ABS ( TEST)-. 000001)13, 13, 55 
10  XM=XM+DXM 
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OXM=*DXM/ 10. 

GO  TO  11 

7  I F ( TEST ) 11«13*  12 

12  XN*XM+DXM 
DXM*DXM/ 10. 

GO  TO  11 

13  RETURN 
END 

FUNCTION  TABLE1(F1,F2,F3,NPTS) 
DIMENSION  Fit  100) ,F2( 100) 
IF(F2( 1 )  —  F  2 ( NPTS ) 1230,230,235 
235  DO  240  K*1,NPTS 
I*K 

I F ( F2  t  I  )-F3 ) 30 , 20 , 240 
240  CONTINUE 
230  DO  10  K= 1 , NPT S 
I  *K 

I F ( F  2  t I ) -F  3 ) 10,20,30 
10  CONTINUE 
20  TABLE 1=F 1 ( I ) 
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GO  TO  AO 
30  IFII-1)1,1.2 

1  TABLE1*F1(I) 

GO  TO  40 

2  A 1  *  F-  2  C  1-1 ) 

A2-F11 I-l) 

3  TABL£1={F1{I)-A2)*(F3-AI)/(F2CI)-A1)+A2 
40  CONTINUE 

RETURN 

END 

END-OF-DATA  ENCOUNTERED  ON  SYSTEM  INPUT  FILE. 
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J«0 
KE  Y  =  0 
DXM= 1 • 

55  I F  t  J-50 ) 5  ,  13,13 

5  FMI  =  SQRT(GPl/GMl)*ATAN(SQRT{GMi/GPl*(XM*XM-l. ) ) )  - AT  AN ( SORT ( XM»XM 

1-1.)) 

TEST^FM I-ANGLE 
IF(KEY)4,4,3 
A  XM-XM-DXM 

IFITESTJ8, 13,9 
9  KEY=  1 
GO  TO  5 
8  KE  Y  =  2 
GO  TO  5 

3  GO  TO  (6,7) , KEY 

6  I F ( TEST ) 10 , 13 , 1 1 
11  XM=XM— OXM 

J  =  J+1 

IF (A8S( TEST) -.000001) 13, 13, 55 
10  XM»XM>DXM 
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OXM=DXM/iO. 

GO  TU  11 

7  IF(TESf) 11,13,12 

12  XM=XM+DXM 
DXM=DXM/10. 

GO  TO  11 

13  RETURN 
END 

FUNCTION  TABLE1(F1,F2»F3,NPTS) 

DIMENSION  F 1 ( 100) , F2(  100) 

IF(F2( 1 )-F2 ( NPTS  > ) 230, 2  30,2 35 
235  DO  240  K=1,NPTS 
1  =K 

IFIF21 I )-F3) 30,20,240 
240  CONTINUE 
230  DO  10  K=  1  *  NPT  S 
I  =K 

1FIF21  I  ) -F  3 )  10,20,30 
10  CONTINUE 
20  TABLE1=F l i I ) 
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GO  TO  AO 
30  IF( 1-1)1,  1,2 

1  T ABLE I  =  F i  {  I  ) 

GO  TO  AO 

2  A 1  *F2  (  1-1  ) 

A2  =  F  1  (  I-l) 

3  TABLE1-(F1{I)-A2)*{F3-Al)/{F2(1)-Ai)+A2 
40  CONTINUE 

RE  FURN 
END 

ENO-OF-OATA  ENCOUNTERED  CN  SYSTEM  INPUT  FILE. 


